function env = calc_env(y, sr, ef_type, flpo, fcut)
% function env = calc_env(y, sr, ef_type, flpo, fcut)
%
% This function implements two simple envelope followers.  The first
% uses the analytic signal (via the Hilbert transform), while the second
% uses a simple low-pass filtering scheme.
%
% INPUT:
%    y       [vec] - Signal to calculate envelope of
%    sr      [num] - Sampling rate
%    ef_type [str] - Type of follower to use ('hilbert' or 'fir')
%    flpo    [num] - Low-pass filter order (if 0, don't filter at all)
%    fcut    [num] - Low-pass filter cutoff frequency (hertz)
%
% OUPUT:
%    env     [vec] - The envelope
%

% Copyright (C) 2010-2030 Graham Grindlay (grindlay@ee.columbia.edu)
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

if nargin < 5
   fcut = 30;
end

if nargin < 4
   flpo = 3;
end

if nargin < 3
   ef_type = 'hilbert';
end

% if only data is passed, default to hilbert follower with no filtering
if nargin < 2
   warning('no sampling rate passed in so defaulting to no filtering');
   ef_type = 'hilbert';
   flpo = 0;
end

% calculate envelope
switch ef_type
    case 'hilbert'
        % reverberant hilbert envelope follower
        yh = hilbert(y);
        
        % square the hilbert transform and take the
        % +ve square root to get the estimated envelope
        env = sqrt(yh.*conj(yh));
        
        % low-pass filter (optional)
        if flpo > 0
            [b,a] = butter(flpo, fcut/(sr/2));
            %envr(c,:) = filter(b, a, envr(c,:));
            env = filtfilt(b, a, env);
        end
        
    case 'fir'
        % simple FIR-lowpass envelope follower
        [b,a] = butter(flpo, fcut/(sr/2));
        
        % reverberant signal
        yl = filter(b, a, 2*y.^2);
        env = sqrt(abs(yl));
        
    otherwise
        error('unknown ef_type!');
end

% clip negative envelope to 0
env(env<0) = 0;

